function [results] = mimo_channel_analysis(varargin)
% MIMO信道容量分析
% 功能：实现不同天线配置的MIMO容量分析、信道条件数分析、自由度分析
% 输入参数：
%   - SNR_dB: 信噪比范围 (dB), 默认 -5:5:30
%   - antenna_configs: 天线配置数组, 默认 [1, 2, 4, 8, 16]
%   - num_realizations: 蒙特卡洛仿真次数, 默认 100
%   - carrier_freq: 载波频率 (Hz), 默认 2.4e9
%   - speed_kmh: 移动速度 (km/h), 默认 50
% 输出：
%   - results: 结构体，包含容量分析结果和可视化数据
%
% 作者：周勇
% 日期：2024年

%% 参数解析
p = inputParser;
addParameter(p, 'SNR_dB', -5:5:30, @isnumeric);
addParameter(p, 'antenna_configs', [1, 2, 4, 8, 16], @isnumeric);
addParameter(p, 'num_realizations', 100, @isnumeric);
addParameter(p, 'carrier_freq', 2.4e9, @isnumeric);
addParameter(p, 'speed_kmh', 50, @isnumeric);
parse(p, varargin{:});

params = p.Results;

%% 添加路径
addpath('../Common');

%% 获取颜色定义
colors = color_definitions();

%% 基本参数计算
SNR_linear = 10.^(params.SNR_dB/10);
speed_ms = params.speed_kmh * 1000 / 3600;
wavelength = 3e8 / params.carrier_freq;
max_doppler = speed_ms / wavelength;

fprintf('=== MIMO信道容量分析 ===\n');
fprintf('SNR范围: %d 到 %d dB\n', min(params.SNR_dB), max(params.SNR_dB));
fprintf('天线配置: %s\n', mat2str(params.antenna_configs));
fprintf('蒙特卡洛仿真次数: %d\n', params.num_realizations);

%% 初始化容量数组
fading_capacity_csi = zeros(length(params.SNR_dB), length(params.antenna_configs));
fading_capacity_no_csi = zeros(length(params.SNR_dB), length(params.antenna_configs));
condition_numbers = cell(length(params.antenna_configs), 1);
effective_dofs = zeros(length(params.antenna_configs), 1);

%% 不同天线配置的容量分析
for ant_config = 1:length(params.antenna_configs)
    Nt = params.antenna_configs(ant_config);
    Nr = params.antenna_configs(ant_config);
    
    fprintf('\n仿真 %dx%d MIMO系统...\n', Nt, Nr);
    
    % 存储条件数样本
    condition_samples = zeros(params.num_realizations, 1);
    
    for snr_idx = 1:length(params.SNR_dB)
        snr = SNR_linear(snr_idx);
        
        % 蒙特卡洛仿真
        capacities_csi = zeros(params.num_realizations, 1);
        capacities_no_csi = zeros(params.num_realizations, 1);
        
        for real = 1:params.num_realizations
            % 生成瑞利衰落信道
            H_fading = sqrt(0.5) * (randn(Nr, Nt) + 1i * randn(Nr, Nt));
            
            % 已知CSI容量 (注水算法)
            [U, S, V] = svd(H_fading);
            singular_values = diag(S);
            capacities_csi(real) = water_filling_capacity(singular_values, snr);
            
            % 未知CSI容量 (等功率分配)
            capacities_no_csi(real) = log2(det(eye(Nr) + (snr/Nt) * H_fading * H_fading'));
            
            % 存储条件数
            if snr_idx == 1 % 只在第一个SNR点计算条件数
                condition_samples(real) = cond(H_fading);
            end
        end
        
        % 遍历容量 (平均值)
        fading_capacity_csi(snr_idx, ant_config) = mean(capacities_csi);
        fading_capacity_no_csi(snr_idx, ant_config) = mean(capacities_no_csi);
    end
    
    % 存储条件数分布
    condition_numbers{ant_config} = condition_samples;
    
    % 计算有效自由度 (基于第一个信道实现)
    H_sample = sqrt(0.5) * (randn(Nr, Nt) + 1i * randn(Nr, Nt));
    [~, S_sample, ~] = svd(H_sample);
    singular_values_sample = diag(S_sample);
    effective_dofs(ant_config) = sum(singular_values_sample.^2) / max(singular_values_sample)^2;
end

%% 可视化结果
figure('Name', 'MIMO信道容量分析', 'Position', [150, 150, 1200, 800]);

% 已知CSI容量
subplot(2,2,1);
for ant_config = 1:length(params.antenna_configs)
    plot(params.SNR_dB, fading_capacity_csi(:, ant_config), ...
         ['-', colors(ant_config)], 'LineWidth', 2);
    hold on;
end
grid on;
xlabel('SNR (dB)');
ylabel('容量 (bps/Hz)');
title('MIMO容量 (已知CSI)');
legend(arrayfun(@(x) sprintf('%dx%d', x, x), params.antenna_configs, 'UniformOutput', false));

% CSI知识的影响 (以4x4为例)
subplot(2,2,2);
ant_idx = find(params.antenna_configs == 4, 1);
if ~isempty(ant_idx)
    plot(params.SNR_dB, fading_capacity_csi(:, ant_idx), 'b-', 'LineWidth', 2);
    hold on;
    plot(params.SNR_dB, fading_capacity_no_csi(:, ant_idx), 'r--', 'LineWidth', 2);
    grid on;
    xlabel('SNR (dB)');
    ylabel('容量 (bps/Hz)');
    title(sprintf('%dx%d MIMO: CSI影响', 4, 4));
    legend('已知CSI', '未知CSI');
end

% 条件数分布
subplot(2,2,3);
ant_idx = 3; % 4x4 MIMO
if ant_idx <= length(params.antenna_configs)
    histogram(condition_numbers{ant_idx}, 20, 'Normalization', 'pdf');
    grid on;
    xlabel('条件数');
    ylabel('概率密度');
    title(sprintf('%dx%d MIMO条件数分布', ...
        params.antenna_configs(ant_idx), params.antenna_configs(ant_idx)));
end

% 有效自由度
subplot(2,2,4);
plot(params.antenna_configs, effective_dofs, 'o-', 'LineWidth', 2);
grid on;
xlabel('天线数量');
ylabel('有效自由度');
title('MIMO有效自由度');

%% 结果打包
results = struct();
results.SNR_dB = params.SNR_dB;
results.antenna_configs = params.antenna_configs;
results.fading_capacity_csi = fading_capacity_csi;
results.fading_capacity_no_csi = fading_capacity_no_csi;
results.condition_numbers = condition_numbers;
results.effective_dofs = effective_dofs;
results.params = params;

fprintf('\nMIMO信道容量分析完成！\n');

end

%% 辅助函数
function capacity = water_filling_capacity(singular_values, snr)
% 注水算法容量计算
num_modes = length(singular_values);
total_power = length(singular_values); % 总功率约束
noise_power = 1/snr;

% 寻找最优水位线
water_level = 0;
for k = 1:num_modes
    temp_level = (total_power + noise_power * sum(1./singular_values(1:k).^2)) / k;
    temp_powers = max(0, temp_level - noise_power./singular_values(1:k).^2);
    
    if all(temp_powers >= 0)
        water_level = temp_level;
    else
        break;
    end
end

% 计算容量
capacity = 0;
for i = 1:num_modes
    if singular_values(i)^2 > 0
        power_alloc = max(0, water_level - noise_power/singular_values(i)^2);
        capacity = capacity + log2(1 + snr * power_alloc * singular_values(i)^2);
    end
end
end